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ABSTRACT 

We have studied the dynamical evolution of rotating globular clusters with direct 
A^-body models. Our initial models are rotating King models; we obtained results 
for both equal-mass systems and systems composed out of two mass components. 
Previous investigations using a Fokker-Planck solver have revealed that rotation has 
a noticeable influence on stellar systems like globular clusters, which evolve by two- 
body relaxation. In particular, it accelerates their dynamical evolution through the 
gravogyro instability. We have validated the occurence of the gravogyro instability with 
direct A^-body models. In the case of systems composed out of two mass components, 
mass segregation takes place, which competes with the rotation in the acceleration of 
the core collapse. The "accelerating" effect of rotation has not been detected in our 
isolated two-mass A^-body models. Last, but not least, we have looked at rotating N- 
body models in a tidal field within the tidal approximation. It turns out that rotation 
increases the escape rate significantly. A difference between retrograde and prograde 
rotating star clusters occurs with respect to the orbit of the star cluster around the 
Galaxy, which is due to the presence of a "third integral" and chaotic scattering, 
respectively. 

Key words: A/'-body simulations - rotating star clusters. 



1 INTRODUCTION 

This paper is devoted to an intriguing problem in the dy- 
namics of globular clusters: How an overall rotation of a 
globular cluster affects its dynamical evolution. There is a 
wealth of observational data which indicates that typical 
globular clusters are rotating. The first method for the de- 
tection of rotation in globular clusters dates back in the first 
half of 20th century. Observations showed that "flattening" 
(i.e. a deviation of spherical symmetry) is a common fea- 
ture of globular clusters. The first ellipticity measurements 
have been done by H. Shapley in the 1920s, see Shapley 
(1930). Subsequent ellipticity measurments of Galactic and 
extragalactic globular clusters have been done, among oth- 
ers, by Geisler & Hodge (1980), Frenk k Fall (1982), Geyer 
et al. (1983), White & Shawl (1987), Staneva, Spassova & 
Golev (1996). The ellipticities of Galactic globular clusters 
are now given in the Harris ca talogueQ cf. Harris (1996). 
A second method for the detection of rotation of globular 
clusters relies on radial velocity measurements of individ- 
ual stars in globulars (Wilson & Coffeen 1954, Mayor et al. 
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1984, Meylan k Mayor 1986, Gebhardt et al. 1994, 1995, 
2000, Reijns et al. 2006 and references therein). The mea- 
sured radial velocities can be plotted against the position 
angle. The amplitude of a sine curve which is fitted through 
the data points serves as a rough measure (i.e. a lower limit) 
for the degree of rotation of the star cluster. Proper motions 
of individual stars have been measured as well and used for 
the detection of rotation in globular clusters (e.g. Peterson 
& Cudworth 1994, van Leeuwen et al. 2000 or Anderson & 
King 2003. The authors of the latter two references conclude 
for uj Cen and 47 Tuc that their axes of rotation are con- 
siderably inclined to the plane of the sky). Table 7.2 in the 
review by Meylan & Heggie (1997) gives already an overview 
of the degrees of rotation of Galactic globular clusters. For 
instance, Gebhardt et al. (1995) estimated Vvot/criD — 0.2 
for the evolved clusters M15 and 47 Tuc. However, globu- 
lar clusters in the Large and Small Magellanic Clouds are 
younger and markedly more elliptical than typical Galactic 
globulars (see the above references), which suggests higher 
values of Vrot/criD- 

From the theoretical side, the first published study of a 
rotating system of A/" = 100 bodies is the work by Aarseth 
(1969) which aimed at the dynamical evolution of clusters 
of galaxies. He already detected a flattening caused by rota- 
tion. Subsequent attempts to study rotating globular clus- 
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ters with direct A/'-body models have been made by Wielen 
(1974) for systems with particle numbers A/" = 50 — 250 and 
Akiyama & Sugimoto (1989) for N = 1000. Wielen con- 
cludes that "a slow overall rotation of a globular cluster 
does not significantly affect the dynamical evolution", and 
Akiyama & Sugimoto do not yet clearly demonstrate numer- 
ically the existence of the gravogyro instability (see below) 
"because of statistical noise due to a small number of par- 
ticles". More recent numerical studies of rotating globular 
clusters have been done using the Fokker-Planck [hereafter: 
FP] method, see the PhD thesis of Einsel (1996), Einsel & 
Spurzem (1999) [hereafter ES99], the follow-up papers by 
Kim et al. (2002, 2004) and the PhD thesis of Kim (2003) 
[hereafter K02, K03, K04]. Under certain approximations, 
the evolution of a phase space distribution function which 
describes the macroscopic state of a globular cluster, is fol- 
lowed by the numerical solution of the FP equation. The ad- 
vantage of such approximate models is their computational 
speed, i.e. one can do large parameter studies in a reason- 
able time. For the most recent FP models including rotation, 
which are compared directly to observations of M5, M15, 47 
Tuc and u Cen, see Fiestas et al. (2006)0 

The present study is again concerned with A^-body 
models of rotating globular clusters. Based on Newton's law 
of gravitation, A^-body modelling is probably the most di- 
rect way to investigate the dynamics of stars in globular 
clusters. However, it is computationally expensive and time- 
consuminglfl On the other hand, direct A^-body modelling 
particularly allows to test the validity of assumptions used 
in statistical modelling of globular clusters, like FP models. 

This paper is organized as follows: Section 2 describes 
our initial models and the numerical method. In Section 3 
we will have a closer look at the effect of rotation on the 
core collapse in isolated equal-mass models. It is well-known 
that the core collapse is caused by the gravothermal insta- 
bility: The negative heat capacity of the self- gravitating core 
of the star cluster leads to an amplification of temperature 
gradients if the halo of the star cluster is extensive enough. 
The core collapse is stopped by the formation of a few bina- 
ries in the core of the star cluster by three-body processes 
(Giersz & Heggie 1994b, 1996). These produce heat at a rate 
slightly higher than the rate of energy loss due to two-body 
relaxation, thereby cooling the cluster core which has nega- 
tive heat capacity. If the star cluster rotates, we numerically 
demonstrate in Section 3 that an additional instability oc- 
curs: The gravogyro instability, which accelerates the core 
collapse in our equal-mass models. In Section 4, which is 
dedicated to tidally limited models, we additionally stud_^ 
the escape of stars from rotating star clusters, i.e. mass losi 



A database of rotating star cluster models exists at the ARI: 
http: / /www. ari.uni-heidelberg.de/clusterdata/ 

Therefore A^-body codes for massively parallel computers have 
been written. Moreover, special-purpose hardware systems like 
GRAPE boards have been in use for A/'-body simulations since 
the early 1990s (see http: //grape. astron. s. u-tokyo. ac. j p/grape/ ). 
A new development is a technology based on FPGAs (Field Pro- 
grammable Gate Arrays). However, such hardware is not used in 
the present study. 

^ The term "mass loss" is not used in the usual sense (due to 
stellar evolution), but in the sense of a mass loss of the whole star 
cluster through escaping stars (i.e. due to dynamical interactions). 



across the tidal boundary. In Section 5, which treats rotat- 
ing systems composed out of two mass components, we study 
the interplay between rotation and mass segregation. 



2 INITIAL MODELS / NUMERICAL METHOD 

We employ generalized King models with rotation. The dis- 
tribution function is given by 



f{E, J.) = C- 



exp 



• exp 



(1) 



which depends on the integrals of motion E (energy) and Jz 
(z-component of angular momentum), where $t is the po- 
tential at the outer boundary of the model (i.e. at the radius 
where the density approaches zero) , ax is the King velocity 
dispersion (which is not the central velocity dispersion) and 
Qo is close to the angular speed in the cluster center (which 
has nearly solid body rotation) . In the distribution function 
([1]), the first two factors on the right-hand side represent 
just a usual King model. The last factor, which depends on 
J2, introduces rotation into the model, i.e. a difference be- 
tween positive and negative Jz ■ It also allows for anisotropy, 
i.e. we have fXr = cr^ / cr^ for the velocity dispersions in 
radial, vertical and tangential directions, respectively. Note 
that models where / oc exp(— 7J2) have been first used by 
Prendergast & Tomer (1970) and Wilson (1975) to model 
elliptical galaxies. However, this term in combination with a 
King model has been first used in Goodman (1983, unpub- 
lished), then by Lagoute & Longaretti (1996), Longaretti & 
Lagoute (1996), Einsel (1996), ES99, K02, K03 and K04, 
Fiestas (2006) and Fiestas et al. (2006). 

The rotating King model based on the distribution func- 
tion ([T|) possesses rotational symmetry and differential ro- 
tation. It has two dimensionless free parameters: The King 
parameter Wq — — ($0 — ^t)/crK, where $0 is the central po- 
tential, and the rotation parameter ujq — y^d/ {AnGpc) • ^0, 
where G is the gravitational constant and pc is the central 
density. It reduces to a King model in the limit cjo = 0. 
ES99 give the ratio A = Trot/^kin of rotational to kinetic 
energy in percent in their Table 1 for rotating King models 
with Wo = 6. The ratio of the root mean squared rota- 
tional velocity to the velocity dispersion is then given by 
Vrot/criD — v^3A/ (1 — A), which may be compared with ob- 
servations. Furthermore, it may be of interest to note that 
between 0.7 < uq < 0.8 there is a transition where the initial 
axisymmetric isolated rotating King models become unsta- 
ble to the formation of a bar (P. Berczik, priv. comm.). 

Our models composed out of two mass components 
are initially uniformly mixed and characterized by two ad- 
ditional dimensionless parameters, the stellar mass ratio 
/i = 7712/7711, where mi and 7712 are the individual masses 
of the light and heavy stars, respectively, and the mass frac- 
tion of heavy stars q — M2/M, where M2 is the total mass 
of the heavy component and M is the total mass of the star 
cluster. 

In direct A/'-body models, the main constraint for a rea- 
sonable computing time is the number of particles N. In our 
simulations, we use N = 5000 for all models with /x ^ 10 
(including all equal-mass models) and N = 32000 for mod- 
els with jjL — 2b and = 50. A typical simulation for one of 
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our equal-mass models took about three days on a 3 GHz 
Pentium PC, until the core collapse time was reached. 

The codes used are nbody6++ (e.g. Spurzem 1999) 
and NBODY6 (see Aarseth 1999, 2003 for an overview). The 
former code is a variant of the latter modified in order to 
run on massively parallel computers. Both codes share, aside 
from parallelization, the same fundamental features and 
yield comparable results. A fourth-order Hermite scheme, 
applied first by Makino & Aarseth (1992), is used for the 
direct integration of the Newtonian equations of motion of 
the A^-body system. The codes use adaptive and individual 
time steps, which are organized in hierarchical block time 
steps, the Ahmad-Cohen neighbor scheme (Ahmad & Cohen 
1973), Kustaanheimo-Stiefel regularization of close encoun- 
ters (Kustaanheimo & Stiefel 1965) and Chain regularization 
(Mikkola & Aarseth 1990, 1993, 1996,1998). 

All models done with nbody6H — h have been calcu- 
lated on single-processor machines at the ARI, except the 
two-mass models TM9-TM12 (see section 5). Some models 
(from series EM2, see Section 3) have been calculated using 
Aarseth's original code nbody6. For some models, we used 
averaging over several runs with the same particle number 
N but different initial values for positions and velocities. 
This was done in order to obtain a better statistical quality, 
cf. Section 2.2 of Giersz & Heggie (1994a) for a discussion 
of this approach. Averaging increases the computing time 
for gravitational forces of a system of N bodies on a linear 
scale, whereas increasing the particle number N increases 
the computing time for gravitational forces oc N"^ in the 
limit of high N. Since we consider in the present work pro- 
cesses acting on the relaxation time scale, an extra power 
of N contributes to the scaUng of the total computing time 
of a model. In our averaged models, the physical quantities 
shown in the plots are determined as the arithmetic mean of 
those values of the considered quantity, which resulted from 
the individual runs. 

The quantities on the vertical axis of our figures are 
shown in dimensionless A/'-body units, i.e. G — M — —AE — 
1, where G is the gravitational constant, M is the total 
mass of the system and E the total energy, see Heggie & 
Mathieu (1986). The resulting length unit is the virial ra- 
dius rv = GM'^ /{-AE), which is of the order of the half- 
mass radius of a globular cluster, see, for instance. Table 
1 in Giirkan et al. (2004). The resulting time unit is then 
given by tv — [GM lr\^~^^'^ . The time on the horizontal 
axis of most figures is shown in units of the initial half-mass 
relaxation time, i.e. the trh means trh(O). For the half-mass 
relaxation time, we adopted the expression 



Model 






N 


Averaging 


tec /trh 


EMla 


3 


0.0 


5K 


no 


10.68 


EMlb 


3 


0.3 


5K 


no 


10.10 


EMlc 


3 


0.6 


5K 


no 


10.35 


EM2a 


6 


0.0 


5K 


4 


7.51 ± 0.52 


EM2b 


6 


0.3 


5K 


4 


7.45 ± 0.25 


EM2c 


6 


0.6 


5K 


4 


7.03 ± 0.19 


EM2d 


6 


0.9 


5K 


3 


6.19 ± 0.35 



Table 1. The initial isolated equal-mass models. For the models 
of series EM2 we also give the run-to-run variation (7n—i in tec /trh 
in the last column. 



sion) and the mean stellar mass, respectively. It may be of 
interest to note that ([2]) can be written in terms of tv as 



tr 



3 • 15.41n(7Ar) 



tv. 



(4) 



For (N = 5000, 7 = 0.11) we have trh ^ 152 ty.The crossing 
time at the virial radius is given by tcr = 2rv/vm — 2^/2 tv- 
For the data evaluation, we employed boxcar smooth- 
ing to reduce noise from our plots. We usually employed 
smoothing widths of 5 — 20 ty. 



3 ISOLATED EQUAL-MASS MODELS 



In this section, we study the effect of rotation on the core col- 
lapse in isolated A^-body models (i.e. without tidal field) and 
find, that these models suggest the occurence of the gravo- 
gyro instability, which has been found earlier in FP models. 
Table [1] shows our initial isolated equal-mass models. 



3.1 Core collapse 

For the theory of the gravothermal instability, which leads 
to the core collapse of star clusters we refer to the origi- 
nal works of Antonov (1962), Lynden-Beh & Wood (1968) 
and Hachisu et al. (1978a/b). It is nowadays also described 
in many textbooks such as Binney & Tremaine (1987) or 
Heggie & Hut (2003). 

Figure [T] shows the time evolution of the core radius 
for the models of series EM2. In this series, one very fast 
rotating model with uoq = 0.9 is included. The core radius is 
defined as 



trh 



8^ [y^2\E\/Mrv]^ N 



(2) 



3 15.4^2 71^2 ln(7Ar)' 

with 7 = 0.11. The above definition is based on Equation 
(5) of the paper by Spitzer & Hart (1971) with 



3A^ 



m 

M 



M 
~N 



(3) 



for the mean value of the particle density inside rv , the root 
mean squared stellar velocity (i.e. the 3D velocity disper- 



rd\'^pi 



(5) 



where fi is the position of the iih star, fd are the coor- 
dinates of the density center of the star cluster and pi is 
the local density at the position of the iih star. The latter 
two quantities are calculated according to the description in 
Casertano & Hut (1985) using the distance to the fifth near- 
est neighbor. It can be seen in Figure [T]that the core radius 
decreases by two orders of magnitude and that the core col- 
lapse occurs earlier for the faster rotating models. Note that 
the four models of series EM2 do not exactly have the same 
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Figure 1. Time evolution of the core radius for the isolated mod- 
els EM2a-d. 



Gravothermal instability 


Gravogyro instability 


drcr'^ / 


drOO / 


Negative specific heat 


Negative specific moment of 


capacity 


inertia 


Energy transport 


Angular momentum transport 


Heat conduction 


Viscosity 



Table 2. A comparison between gravothermal and gravogyro in- 
stability. 



initial core radius in the beginning, i.e. the core radii differ 
by factors up to approximately two. This indicates that the 
models are indeed not in all respects exactly comparable, 
cf. Figures 14/15 of K02. They have slightly different den- 
sity distributions corresponding to their degree of flattening. 
However, since the core shrinks by two orders of magnitude 
during the core collapse, a factor of two in the difference of 
the initial core radii seems to be negligible. 

Figure [2] shows the time evolution of the Lagrangian 
radii for the models of series EM2. The inner Lagrangian 
radii shrink and the outer ones expand. 

The last column of Table [T] contains the core collapse 
times for the isolated equal-mass models. The core collapse 
time is defined as the time when the core radius Vc reached 
its first sharp minimum. For the averaged runs the given 
time is the arithmetic mean of the core collapse times which 
we determined for each single run. We also give the standard 
deviations dn-i as a measure of the run-to-run variation of 
the core collapse times. 



3.2 Gravogyro instability revisited - I 

There is obviously a need to explain the accelerated core 
collapse of the rotating models in Section 3.1. A theoreti- 
cal model, which explains such an effect, has been proposed 
as early as in the 1970s by Inagaki & Hachisu (1978) and 
Hachisu (1979, 1982). They called it the gravogyro insta- 
bility. This instability is expected in self- gravitating and 
rotating systems and has been derived for self- gravitating 
cylinders of infinite length in z-direction: If we equate the 




t/U 

Figure 3. Time evolution of the z-component of the specific an- 
gular momentum of the particles within the Lagrangian radii for 
the fast rotating isolated model EM2c. The isolated points are 
artefacts and stem (like the interruptions in the 1% curve) from 
the sign change of jz • 



centrifugal force to the gravity in the equatorial plane (z = 
■ -\- y^) and neglect the pressure gradient, we have 



2GMr 
R 



(6) 



where Mr is the mass per unit length in 2;-direction con- 
tained within the radius R in cylindrical coordinates and jz 
is the z-component of specific angular momentum. We read- 
ily obtain Sjz/jz — SR/R. Using the equation jz = R^ou in 
addition, we obtain Sjz — u;R6R and 6jz = 2ooRSR-\- R^Soo. 
The combination of these relations yields the following lin- 
ear differential relation between specific angular momentum 
jz and angular speed u: 



(7) 



A negative specific moment of inertia occurs in this relation, 
cL Inagaki & Hachisu (1978). 

If we let an initial model evolve, which exhibits a radial 
gradient of the angular speed and has reasonable rotation 
curve, angular momentum is transported outwards on the 
relaxation time scale by the stellar dynamical analog of vis- 
cosity. The core of the star cluster contracts because of a 
deficit in the centrifugal force. The angular speed in the 
core increases according to relation (jT]). A runaway depar- 
ture from the initial state takes place. This cycle is what 
they denoted as the gravogyro instability. For a comparison 
between gravothermal and gravogyro instabilities see Table 

n 

One may wonder: Does the gravogyro catastrophe re- 
ally occur in isolated A/'-body models? Figure [3] shows the 
time evolution of the ^-component of specific angular mo- 
mentum within the Lagrangian radii for the fast rotating 
isolated model EM2c. The z-component of specific angular 
momentum of each star is calculated according to 



N-body Models of Rotating Globular Clusters 5 




Figure 2. Time evolution of the Lagrangian radii, which contain 1%, 2%, 5%, 10%, 20%, 30%, 40%, 50%, 70% and 90% of the initial 
total mass of the star cluster for the isolated equal-mass models EM2a-d. The inner Lagrangian radii contract, while the outer ones 
expand in the core collapse phase. 
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jz [{r^ - Td) X iJ^]^, (8) 

where = is the position of a particle, = 

{v^x,v^y,v^z) is its velocity and fd = (xd,yd,Zd) is the den- 
sity center of the star cluster. It is then summed over all 
particles within a given Lagrangian radius and divided by 
the number of particles inside that radius. One notes that 
angular momentum does indeed diffuse from the inner parts 
of the cluster to the outer parts, as time proceeds. 
Now we study the evolution of the angular speed oo. The 
gravogyro instability implies that jz goes down while u in- 
creases. Figure |4] shows the time evolution of the average 
angular speed of the particles within the Lagrangian radii. 
The angular speed of one particle is calculated according to 

M ' 

with the same notations as in equation (|8]), but i?^ = 
a/ {Xi, — XdY + (y* — VdY is the radius with respect to the 
density center in cylindrical coordinates. It is then summed 
over all particles inside a given Lagrangian radius and di- 
vided by the number of particles within that radius. One can 
observe an increasing average angular speed inside those La- 
grangian radii, which show a decrease in the z-component of 
angular momentum. Thus, it is shown that gravogyro effects 
appear in our isolated models. 



4 TIDALLY LIMITED EQUAL-MASS MODELS 

As a step towards more realistic models of globular clusters, 
we investigate in this section the effects of the tidal field of 
the Galaxy on the dynamical evolution of rotating globular 
clusters. The implementation of the tidal field of the Milky 
Way Galaxy within the tidal approximation used in nbodyG 
and nbodyGH — h is in detail described in Appendix A. Dif- 
ferent escape criteria are discussed in Appendix B. In our 



Model 






N 


Averaging 


^cc/^rh 


EM3a 


6 


0.0 


5K 


4 


7.11 ± 0.37 


EM3b 


6 


0.3 


5K 


4 


6.37 ± 0.36 


EM3c 


6 


0.6 


5K 


4 


5.42 ± 0.33 


EM3d 


6 


0.9 


5K 


4 


3.97 ± 0.40 


EM4b 


6 


-0.3 


5K 


4 


6.50 ± 0.33 


EM4c 


6 


-0.6 


5K 


4 


5.62 ± 0.07 


EM4d 


6 


-0.9 


5K 


4 


4.88 ± 0.28 



Table 3. The initial tidally limited equal-mass models. The rotat- 
ing models of series EM3 are prograde rotating while the models 
of series EM4 are retrograde rotating with respect to the orbit of 
the star cluster around the Galaxy. We also give the run-to-run 
variation (7n—i in tccl'^rh i^ the last column. 

models, the star cluster initially completely fills its Roche 
lobe in the tidal field of the Galaxy. In other words, the 
density of the initial rotating King model approaches zero 
at the physical tidal cutoff radius. The tidal radius is defined 
as the distance of the star cluster center to the Lagrangian 
points L1/L2. 

Table [3] shows the initial tidally limited equal-mass 
models. The model of series EMS are prograde rotating, i.e. 
the stars move around the z-axis of the cluster-centered co- 
ordinate system in the same sense as the star cluster moves 
around the Galaxy. The models of series EM5 are retrograde 
rotating, which we denoted by a negative sign of the rota- 
tion parameter cjo for clarity. The only difference between 
the rotating initial models of series EMS and EM4 is that 
we reversed the sign of all initial velocity vectors. 

4.1 Core collapse 

Figure [6] shows the time evolution of the core radius for the 
prograde rotating tidally limited models EMSa-d. A com- 
parison of Figures [1] and [6] and the core collapse times in 
Tables [1] and [3] shows, that the tidal field further accelerates 
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solid: N-body, N=5000, radius cutoff (A, Ernst) 

lower lines: prograde, upper lines: retrograde 
dotted; FOPAX, energy cutoff (J, Fiestas) 
dashed: N-body, N"=5000, energy cutoff (E. Kim) 
dot-dashed; K-body, N= 16000, energy cutoff (D, Heggie) 








10 

t/t 



15 



rh 



Figure 5. Time evolution of the total mass of the star clusters for the prograde rotating tidally limited models EM3a-d (lower solid 
lines), the retrograde rotating tidally limited models EM4b-c (upper solid lines), FP models with energy cutoff (dotted), two A^-body 
models with N = 5000 and energy cutoff (dashed) and one non-rotating A^'-body model with = 16000 and energy cutoff (dot-dashed). 
The core collapse times of our A^-body models are marked with a cross. 



the core collapse by a significant amount. The reason for this 
behaviour is that rapid mass loss across the tidal boundary 
provides an additional mechanism for the transport and re- 
moval of energy and angular momentum, thus enhancing the 
effects of the gravothermal and gravogyro instabilities. An 
inspection of Table [3] reveals that there is a difference be- 
tween the core collapse times between the models of series 
EMS and EM4, i.e. the prograde rotating models of series 
EMS collapse faster than the retrograde rotating model of 
series EM4. By considering the run-to-run variation an-i in 
the last column of table [S] we see that the statistical sig- 
nificance of this effect is strongest for the two models with 
(jJo — 0.9 as compared with the case of the slower rotating 
models.. Therefore this effect might be related to the forma- 
tion of a bar in the fastest rotating models. On the other 
hand, the difference of the core collapse times between the 
prograde and retrograde rotating models may also be related 
to the difference in the rates of mass loss (see next Section): 
The loss of energy and angular momentum through escaping 
stars is more rapid for the prograde rotating models and thus 



the effects of the gravothermal and gravogyro instabilities, 
which accelerate the core collapse, are stronger. 

4.2 Mass loss 

Figure [5] shows the time evolution of the total mass of the 
globular clusters, where the initial mass is normalized to 
one, for several models, including simulations with the 2D 
FP code FOPAX (with energy cutoff, cf. Fiestas 2006), one 
rotating and one non-rotating A/'-body model by E. Kim 
(with energy cutoff, without the tidal approximation, see 
KOS) and one non-rotating A^-body model by D. C. Heggie 
(with energy cutoff, without the tidal approximation). The 
solid lines are our nbody6H — h results (with radius cutoff 
and the tidal approximation). The different "cutoffs" are 
explained in Appendix B. The initial models are the same 
rotating (or non-rotating) King models. Note that for the 
rotating models, there are two solid lines for each color. The 
lower of these lines corresponds to a prograde rotating model 
(i.e., from series EMS). The upper of these lines corresponds 
to a retrograde rotating model (i.e., from series EM4). The 
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Figure 6. Time evolution of the core radius for the prograde 
rotating tidally hmited models EM3a-d 



escape rate for retrograde rotating models is considerably 
lower, since many stars in retrograde orbits are subject to 
a "third integral" of motion, which restricts their accessible 
phase space and hinders them from escaping, see Fukushige 
& Heggie (2000) and Appendix B. 

Figure [5] also shows, that the faster rotating models suf- 
fer stronger mass loss than the slowly or non-rotating mod- 
els, a fact which has already been shown in the FP models 
of ES99. 

The differences between the models of the same rota- 
tion parameter are due to different implementations of the 
tidal field in the various codes (cf. Appendix B) or due to 
the different modelling techniques in general. For instance, 
the difference between the A/'-body and FP models is the 
tidal approximation used in nbodyG and nbody6+ + , i.e. 
the modified equations of motion with a linear approxima- 
tion of the tidal forces (see Appendix A), a feature, which 
cannot be implemented in an FP code, which relies only 
only on a tidal cutoff. The geometry of the almond-shaped 
tidal boundary within the tidal approximation (see Figure 
IBip differs from the geometry of our star cluster models, 
which are axisymmetric. This may cause another difference 
between A/'-body and FP models. In addition, the difference 
between prograde and retrograde rotating models cannot be 
seen in FP models. Note that Kim et al. (2006) use for their 
comparison A/'-body models with an artificially implemented 
tidal cutoff and no modification of the equations of motion. 

Last but not least, a further inspection of Figure [5] re- 
veals, that core collapse causes an increase in the escape rate, 
i.e. short before the core collapse time the curves steepen 
slightly. The moment of core collapse is marked with a cross 
on each curve corresponding to one of our A^-body mod- 
els. However, our main aim was to show the difference of 
prograde and retrograde rotating models and the fact, that 
rotation significantly increases the rate of mass loss. 

4.3 Gravogyro instability revisited - II 

In this subsection, we briefly look at the gravogyro insta- 
bility in the tidally limited A/'-body models. The mass loss 
through the tidal boundary provides, as noted before, an 
additional mechanism for the transport and removal of an- 




Figure 7. Time evolution of the Lagrangian radii for a single run 
of the fast rotating tidally limited model EM3c. For explanations 
see the text. 




Figure 8. Time evolution of the z-component of the specific an- 
gular momentum of the particles within the Lagrangian radii for 
a single run of the fast rotating tidally limited model EM3c. For 
explanations see the text. 



gular momentum (besides viscosity effects) . Figure [71 shows 
the time evolution of the Lagrangian radii for one run of 
the tidally limited model EM3c. The Lagrangian radii are 
defined with respect to the initial particle number. There- 
fore the particle number within the Lagrangian radii remains 
constant while the system evolves, whereas after some time, 
the outer Lagrangian shells are completely lost due to es- 
caping stars, i.e. the curves for the outer Lagrangian radii 
bend upwards until they reach twice the tidal radius, which 
is the cutoff radius in the adopted escape criterion. Figure [8] 
shows the time evolution of the z- component of specific an- 
gular momentum within the Lagrangian radii for the same 
single run of the tidally limited model EM3c. The curves for 
jz decrease due to escaping stars carrying away angular mo- 
mentum. The curves stop at the point at which the adjacent 
outer Lagrangian shell is lost, such that we always have a 
constant particle number within the Lagrangian radii. Fig- 
ure [9] shows the time evolution of the average angular speed 
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Figure 9. Time evolution of the average angular speed of the particles within the Lagrangian radii for a single run of the fast rotating 
tidally limited model EM3c with Wq =6, uq =0.6. For explanations see the text. 



of the stars within the Lagrangian radii for the same single 
run of the tidally limited model EM3c. After an initial de- 
crease of uj due to mass loss, which is most prominent in the 
outer Lagrangian spheres but reaches deep into the core, oo 
increases more than in the isolated model and up to higher 
radii. A similar effect can be seen in FP models, see Einsel 
(1996), Figures 40 and 42. 



5 ISOLATED TWO-MASS MODELS 

Realistic globular clusters have a mass spectrum. However, 
the main effects of a mass spectrum on the dynamics of 
globular clusters can already be seen in idealized models 
composed out of only two mass components. We study such 
models in this section. Table |4] summarizes our initial models 
with two mass components. All these models are done with 
NBODY6++. We would like to stress that they are isolated, 
and no averaging was used. The total mass fraction of heavy 
stars q was fixed to the value q = 0.1 as for the largest part 
of the parameter space covered in Khalisi (2002) and Khalisi 
et al. (2006). As he notes, this choice was inspired by a result 
of Inagaki & Wiyanto (1984) that a cluster evolves fastest 
for q — 0.1. Following the work of Khalisi, we varied the 
stellar mass ratio /x. The models of series TM9 and TMIO 
with high particle numbers [N—32K) were computed (using 
32-128 processors) at the IBM Regatta p690+ supercom- 
puter "Jump" at the Research Center Jiilich, consisting of 
a total of 1312 Power4+ processors running at a frequency 
of 1.7 GHz, whereas the models of series TMll and TM12 
(also A/'=32K) were computed with the PC Beowulf cluster 
"Hydra" of the ARI, consisting of 10 Dual P4 with 2.2 GHz. 
Each of the core collapse times given in Table |4] is the arith- 
metic mean of three times which were determined for each 
run: When the core radius had its first sharp minimum, the 
minimum potential in the star cluster reached its first sharp 
minimum and the maximum density in the core reached its 
first sharp maximum. 

Figure [10] shows the time evolution of the Lagrangian 
radii for the model TM7a, computed seperately for the 
light /heavy component. While the orbits of the light stars 
expand, the heavy stars sink to the center. As can be seen, 
only the heavy masses are involved in the core collapse. For 
the non-rotating case with q — 0.1, the transition between 
the Spitzer-stable and Spitzer-unstable regimes lies indeed 
between /i = 1.25 and /i = 2, as we checked in our data. 
For details on the Spitzer instability, see the original work 
of Spitzer (1969) or the review in Khalisi (2002). The em- 




Figure 10. Time evolution of the Lagrangian radii of the model 
TM7a, computed seperately for the light masses, neglecting the 
heavy masses and for the heavy masses, neglecting the light 
masses. In particular, this Figure shows the radii, in which 1%, 
90% of the light/heavy masses are contained. 



pirical criterion by Watters et al. (2000), which is based on 
Monte- Carlo simulations, also yields a transition occuring 
between the above values of /i. 

Figure [11] shows the evolution of the maximal density 
in the core for the models of series TM5 - TM8 and TMll 
- TM12. What can immediately be seen in Figure [TTI is the 
acceleration of the core collapse due to mass segregation. 
The higher the stellar mass ratio /i, the shorter is the core 
collapse time. 

Figure [12] shows the acceleration factor of the core col- 
lapse due to mass segregation. For this plot we used the 
equal-mass core collapse times tcc,i of the non-rotating iso- 
lated models EMla and EM2a (see Section 3) and the core 
collapse times tcc,/x from the non-rotating models "a" of se- 
ries TM1-TM12. Khalisi and Spurzem (2001) confirmed the 
relation 



(10) 



which is derived from the relevant time scales and also shown 
in Figure 1121 In the limit of small /i the curve approaches 
the Spitzer-stable regime while in the limit of large /i N- 
dependant effects play a role, which can be seen in more 
detail in the study by Khalisi et al. (2006) (their Figure 6). 
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Figure 11. Time evolution of the maximal density in the core for the models of series TM5 - TM8 and TMll - TM12 
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Table 4. The initial models with two mass components. 
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Figure 12. Acceleration factor of the core collapse due to mass 
segregation. The positions of the crosses (N=5000) and triangles 
(N=32000) are calculated from values in Table 0] for the non- 
rotating models and normalized to values in Table ^for the equal- 
mass models (/i = 1). 



However, in Figure [TT] the acceleration of the core col- 
lapse due to rotation is not detectable. Rather, by also tak- 
ing into account the core collapse times in columns 7 and 
14 of Table IH it seems, that in most models, there is an 
opposite effect, i.e. the rotating models have a delayed core 
collapse. 

Figure [13] shows the time evolution of the Lagrangian 
radii for the models TM7a and TM7c, computed only for 
the heavy masses. A comparison of the curves for the models 
with different uo in Figure [T3l reveals that the inner (< 50%) 
Lagrangian radii for the model (cjq = 0.6) decrease some- 
what slower than those of the non- rotating model {ooo = 0.0). 
This points towards a slowdown of the mass segregation pro- 
cess induced by rotation. 



0,10 r 




Figure 13. Time evolution of the Lagrangian radii of the models 
TM7a and TM7c, computed only for the heavy masses, neglecting 
the light masses. 



Figure [14] illustrates the time evolution of the z- 
component of angular momentum contained within the core, 
the half- mass radius and the 90% Lagrangian radius for the 
two mass species seperately for the isolated model TM6c. For 
each color, the upper curve corresponds to the light masses, 
while the lower curve corresponds to the heavy masses. Note 
that the z-component of the total angular momentum re- 
mains approximately constant, since the model is isolated. 
The core loses angular momentum to the outer parts of the 
cluster. Kim et al. point out, that angular momentum is 
transferred from the high masses to the low masses. We 
cannot prove this trend in Figure [14] As we consider the 
angular momentum contained in the core, it can be seen 
that the high masses and the low masses as well lose angu- 
lar momentum. The angular momentum loss of the heavy 
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Figure 14. Time evolution of the z-component of angular mo- 
mentum contained within the core, the half-mass radius and the 
90% Lagrangian radius, separately for the two mass components 
(for each linestyle: upper curve: light, lower curve: heavy) for the 
model TM6c. The interruptions in the black curves stem from the 
sign change of the plotted quantity. 

masses is a bit slower than that of the light masses which is 
still in accordance with K04 (see their Figure 11). 

At last, we would like to mention that we checked the 
pairing of the dynamically formed binaries in our two-mass 
models. In the Spitzer-stable regime (/x = 1.25) all pair- 
ings occur: Heavy-heavy, light-light and heavy-light. The 
pairings change from time to time. In the Spitzer-unstable 
regime of our parameter space (/x = 2 — 50), we noticed (with 
a few exceptions) only heavy-heavy binaries. 



6 DISCUSSION 

The main aim of this paper was to investigate the effect of 
an overall (differential) rotation on the dynamical evolution 
of globular clusters with direct A^-body models. We would 
like to briefly summarize our results first and then tackle a 
comparison with results obtained by others. 

6.1 Summary 

The difference between rotating and non-rotating equal- 
mass A^-body models of globular clusters is that an overall 
(differential) rotation accelerates their dynamical evolution. 
For instance, the core collapse of equal-mass models is ac- 
celerated for the cases with rotation as compared with the 
non-rotating models. However, this acceleration of the dy- 
namical evolution induced by rotation cannot be seen in our 
isolated two-mass A^-body models. 

As can be seen see in Table 1, the core collapse times 
are larger for the models with King parameter Wq — 3 than 
for those with Wq — 6, which are more concentrated. This is 
a well-known fact, which is noted, for instance, in Quinlan 
(1996) or Giirkan et al. (2004). For the models EMla-c with 
Wq — 3, the acceleration of the core collapse due to rotation 
is not visible within the estimated error margins of the mea- 
surement. Therefore, runs with higher particle numbers are 
required, or, an averaging over several runs with different 



initial configurations of positions and velocities, which has 
been done for the models EM2a-d with Wq = 6 in order to 
damp statistical fluctuations. 

In the tidally limited models, we observe rapid mass loss 
across the tidal boundary, which is stronger for the faster 
rotating models. A tidal field further accelerates the core 
collapse, since it provides an additional mechanism for the 
transport and removal of energy and angular momentum 
through escaping stars, thus enhancing the effects of the 
gravothermal and gravogyro instabilities. 

In the case of tidally limited A'-body models within the 
tidal approximation (see Appendix A), the sense of rota- 
tion plays an imporant role in the escape process. Models, 
where most stars are in retrograde orbits (as compared with 
the sense of rotation of the star cluster around the galaxy), 
have a significantly lower escape rate due to the presence of 
a "third integral", cf. Fukushige & Heggie (2000) and Ap- 
pendix B. 

The gravogyro instability predicted by Inagaki & 
Hachisu (1978) and Hachisu (1979, 1982), as it becomes 
manifest in a decrease of the component of angular mo- 
mentum combined with an increase in the angular speed, 
was found in both isolated and tidally limited equal-mass 
A'-body models. 

In the models with two mass components, mass segre- 
gation takes place, which again accelerates the dynamical 
evolution of the globular cluster and results in a faster core 
collapse, the higher the stellar mass ratio /i is. Our results 
are in fair agreement with results in Khalisi et al. (2006). 

However, if both mass segregation and rotation com- 
pete in the acceleration of the dynamical evolution of the 
globular cluster, the trend that rotation accelerates the dy- 
namical evolution of the stellar system, is no more visible in 
isolated models. Rather, in most models, there seems to be 
an opposite effect: The faster rotating isolated models have 
a delayed core collapse as compared with the non-rotating 
models. A possible explanation is that the rotation slows 
down the mass segregation process. However, the statistical 
quality of our results is limited. 

6.2 Comparisons with other A^-body models 

A few A'-body models of rotating globular clusters have been 
presented in chapter 4 of K03. This chapter is a comparative 
study between A'-body simulations and FP model results; 
see also the recent work by Kim et al. (2006) [hereafter: 
K06]. Kim's initial models are the same rotating King mod- 
els as they are described in Section 2. He used a version 
of Aarseth's original code nbody6 modified "to mimic the 
tidal environment of the clusters modeled with the 2D FP 
equation." Namely, he applied a tidal energy cutoff to his A'- 
body models (see Appendix B). With the same definition of 
relaxation time, the half-mass times of his and our A'-body 
models differ significantly, as can be seen in Figure [5] The 
same is valid for a comparison of the half-mass times of our 
non-rotating A^-body model and the corresponding A^-body 
model by Heggie (also with energy cutoff). The reason for 
these rather large differences is the different treatment of the 
tidal field. The core collapse times from our A^-body models 
agree well with those of Kim's A'-body models. Our results 
are also consistent with A'-body simulations presented in 
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Ardi et al. (2005, 2006) if the same definition of the half- 
mass relaxation time is used. 



6.3 Comparison with FP models 

For a comparative study of A/'-body models of rotating star 
clusters with FP models we refer to chapter 4 of K03 and the 
recent work K06. While the time step in A/'-body simulations 
is a fraction of the orbital time, it is a fraction of relaxation 
time in the FP codes used in ES99 and K02, K03, K04 and 
K06. Therefore, the proper definition of relaxation time is 
crucial for a comparison of A/'-body and FP models. With 
equation Q as definition of the half-mass relaxation time, 
the core collapse times in our equal-mass A^-body models 
(i.e. the models of series EM1-EM4) are shorter than the 
core collapse times in the tidally limited FP models in ES99 
(see their Figure 2) and the isolated FP models in K02 (see 
their Figure 3). A similar trend has been found in K03 and 
K06. 

The different implementations of a tidal boundary in 
A^-body and FP models and the absence of tidal forces in 
FP models have a significant influence on the escape rates 
(see Section 4.2 and Figure [5]). 

In the mdels with two mass components, a difference 
to FP models of K04 occurs: In our isolated models, we do 
not see the effect of acceleration of the core collapse due to 
rotation any more. 



6.4 Future work 

Several questions remained unanswered - the endeavor to 
answer some of the following questions may be the start- 
ing point for future investigations: Does the suspected slow- 
down of the mass segregation process induced by rotation 
also occur in tidally limited models with two mass compo- 
nents? How evolves the average angular speed of the different 
mass components, taken seperately? We need a more inten- 
sive A/'-body study of tidally limited rotating systems with 
a mass spectrum and a better statistical quality, i.e. higher 
particle numbers or extensive averaging. How do stellar evo- 
lution, disk shocking, primordial binaries or a central black 
hole influence the dynamical evolution of rotating globular 
clusters? 
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Figure Al. Equipotential lines in the z=0 plane for the effective 
potential (|A1)) in the case k"^ ~ 1.8 o^q, ujq = G = 1. For the star 
cluster we assumed a Plummer potential, for which we defined a 
"concentration" c = \ogiQ(rt/rpi) = 1 and a total mass M ~ 2.2. 
In these units, we have rt = 1 for the tidal radius and Cl — —3.3 
for the critical Jacobi integral. 

APPENDIX A: THE TIDAL APPROXIMATION 

The mean eccentricity of the known 3D globular cluster or- 
bits in the halo of the Milky Way is rather high, see, e.g., 
Dinescu et al. (1999) or Allen et al. (2006) for a more recent 
work. Therefore, a realistic tidal field acting on a globular 
cluster in the Galactic halo should be time-dependant. How- 
ever, in nbodyG and nbodyGH — h, the star cluster is assumed 
to move around the Galactic center on a circular orbit, which 
makes it possible to use the epicyclic approximation, cf. 
Binney & Tremaine (1987). The corresponding tidal field is 
steady. Its implementation in the A^-body code is as follows: 
As in the restricted three-body problem, one applies a coor- 
dinate transformation to a rotating coordinate system with 
the angular velocity coo, in which both the star cluster cen- 
ter and the Galactic center (i.e., the primaries) are at rest. 
Its origin is the star cluster center, sitting in the minimum 
of the effective Galactic potential. The x-axis points away 
from the Galactic center. The y-axis points in the direction 
of the rotation of the star cluster around the Galactic center. 
The z-axis lies perpendicular to the orbital plane and points 
towards the Galactic North pole. Through this transforma- 
tion into a frame of reference, in which the tidal potential 
is static, centrifugal and Coriolis forces appear according to 
classical mechanics. In addition, tidal terms enter the mod- 
ified equations of motion in the star cluster region. These 
can be derived from an effective potential in the epicyclic 
approximation: 

$eff(x, y, z) = $ci(x, y, z) + i/i^x^ + 

+ 0{xz'^) + const (Al) 

where the coordinates are relative to the star cluster center 
and $ci is the star cluster potential. Since the expansion is 
about the minimum of $eff, where the first partial deriva- 
tives vanish, there are no first-degree terms in the Taylor ex- 
pansion. The term oc xz vanishes because $eff is symmetric 
in z. The term oc arises from a combination of centrifugal 
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Isolated system 


Tidally limited system 


1. E^>0 


1. a > Cl 


2. Ir^ -fd\> rcrit =20 -rv 


2. |r; - rd\>2-rt 



Table Bl. Escape criteria in nbody6 and nbody6++. 

and tidal forces. For an illustration using a Plummer poten- 
tial for the star cluster, see Figure IA1[ which also shows 
the location of the effective potential's nearest equilibrium 
points with respect to the star cluster center, which are cov- 
ered by the approximation: The Lagrangian points Li and 
L2, where V$eff = 0. They are saddle points of the effective 
potential, i.e. the Hessian is neither positive nor negative 
definite. In terms of Oort's constants A and B, we have in 
the solar neighborhood 



0^0 
2 



(0,0, cjo) = (0,0,A-B), 
-4B{A-B), 

- 4cjo = -4A{A - B) < 0, 
47TGpg + 2{A^ - B^) > 0, 



(A2) 
(A3) 
(A4) 
(A5) 



where k and u are the epicyclic and vertical frequencies, 
respectively. The ratio a^^/cJq — 1-8 calculated from Oort's 
constants depends in the general case on the density pro- 
file of the Galaxy. The vertical frequency u can be derived 
from the Poisson equation for an axisymmetric system, see 
Oort (1965), and pg is the local Galactic density, which con- 
tributes to the dominant first term in (lA5l)FI The equations 
of motion in the rotating reference frame are then 



2{ooo X x). 



(A6) 



where the last term on the right-hand side represents the 
Coriolis forces, which cannot be derived as the usual gradient 
of a potential, since they are velocity-dependant. After a 
little bit of vector analysis, the equations of motion read 



f:,+2{A- B)y + 4A{A- B)a 
fy - 2{A - B)x 

- [47rGpg + 2{A^ - B^)] z, 



(A7) 
(A8) 
(A9) 



where (fx, fy, fz) = — V$ci is the force vector from the other 
cluster member stars. 



APPENDIX B: ESCAPE CRITERIA 

In nbodyG and nbodyGH — h, the escape criteria shown in 
Table IBll are implemented. is the energy of a star, 
its Jacobi integral, its position vector, fd the position of 
the star cluster's density center and rv is the virial radius. 
In both the isolated and tidally limited cases, there are a 
necessary and a sufficient criterion for escape implemented 
in nbodyG and nbodyGH — h- Since the sufficient criterion is 
always related to a critical radius, we refer to these escape 
criteria as radius cutoffs. 

^ Actually, such a treatment is only valid in the solar neighbor- 
hood, i.e. in the Galactic disk. 




Figure Bl. The critical equipotential surface at $eff = Cl with 
the same conditions as in Figure lAll We additionally assumed 



The z-axis is pointing upwards. 



The tidal radius is given by 



GM 



4A{A-B) 



1/3 



M 



1/3 



Ro 



(Bl) 



according to King (1962), where G is the gravitational con- 
stant, M is the total mass of the star cluster, A and B are 
Oort's constants Mg is the mass of the Galaxy and Ro is 
the Galactocentric radius of the star cluster's orbit around 
the Galactic center. The critical Jacobi integral is given by 
the value of the effective potential ()A1|) at the Lagrangian 
points Li and L2, 



Cl 



3GM 
2rt 



[G^M^4A{A - B)] 



1/3 



(B2) 



according to Wielen (1972). It is conserved in the rotating 
reference frame, which is centered on the star cluster cen- 
ter. For Ci, > Cl, the equipotential surfaces in Figure Al 
are open and form channels of escape, through which orbits 
can leak out. The last closed surface is almond-shaped (see 
Figure [BT]) . 

Takahashi & Portegies Zwart (1998) note that "In an 
isotropic Fokker-Planck model, one has no choice but to use 
the energy as a criterion for escape." What they refer to is 
the so-called energy cutoff in contrast to the radius cutoff 
of the A/'-body codes. Stars with 



E> Et = $(rt) = 



GM 

rt 



(B3) 



are removed from the system, where Et is the tidal en- 
ergy, i.e. the potential at the tidal radius. This criterion has 
already been used by Chernoff & Weinberg (1990). Stars, 
whose energy exceeds the tidal energy by a small amount, 
are able to cross the tidal radius (if they are on radial orbits). 
In the energy cutoff picture within the FP approximation, 
the lifetime of a star cluster scales exactly with relaxation 
time. However, such an approach may simplify the physics: 
Once the star's energy is higher than the tidal energy it is 
called a "potential escaper" : It has not yet escaped but still 
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C = Cl 




-1 -0.5 0.5 1 



Figure B2. Poincare section of orbits at the critical Jacobi in- 
tegral Cl for the equations of motion (|A7|) - (IA8|) with the same 
conditions as in Figure lAl] at the moment of crossing y = with 
y > 0. The full range of initial conditions is covered. 



cluster around the galaxy, see Fukushige & Heggie (2000). 
The "third integral" restricts their accessible phase space 
and hinders them from escaping, even if the stars have been 
scattered above the critical Jacobi constant. On the other 
hand, the particles on orbits corresponding to the chaotic 
domains in the surface of section bounce back and forth for 
a certain time in a bounded area, called the "scattering re- 
gion" , until they escape through one of the exits, which open 
up around the Lagrangian points Li and L2, for stars with 
a Jacobi constant, which is higher than the critical value. 
For an interesting overview of the physics of chaotic scatter- 
ing see the papers by Aguirre et al. (2001) and Aguirre & 
Sanjuan (2003). 

This paper has been typeset from a TgX/ file prepared 

by the author. 



needs a time of the order of the crossing time to reach the 
tidal radius. Immediate removal of stars which fullfill the 
energy criterion is only reasonable if the orbital timescale is 
negligible compared with the relaxation time. An improved 
approach, which take this dependance on the crossing time 
scale into account, can be found in the paper by Lee & Os- 
triker (1987). 

In an anisotropic FP code, it is possible to use the more 
realistic apocentre criterion by Takahashi et al. (1997) and 
Takahashi & Portegies Zwart (1998): Stars are removed from 
the system, if their apocenter distance Va calculated accord- 
ing to 

J2 = 2rl [E - $(r.)] , (B4) 

exceeds the tidal radius. 

Last but not least, in the case of a radius cutoff and 
the tidal approximation, the picture is similar to the one 
sketched above: Once a star's Jacobi integral has ex- 
ceeded the critical value Cl slightly (i.e. once the star is 
a "potential escaper"), it still needs a certain time to find 
an opening to one of the escape channels in the equipoten- 
tial surfaces shown in Figure Al for the effective potential 
()Aip . This can take many crossing times depending mainly 
on the excess energy, i.e. we have te oc Cl{C^ —Cl)~^ from 
an upper limit on the flux of phase space volume through 
L1/L2 according to MacKay (1990) and Fukushige & Heggie 
(2000). In the radius cutoff picture within the tidal approx- 
imation, the scaling of the lifetime is very subtle and the 
half-mass time scales as tmh oc t^^^, see Baumgardt (2001). 

By integrating the equations of motion ()A7|) - ()A9|) nu- 
merically for orbits at a given Jacobi constant, one can ob- 
tain a Poincare section as shown in Figure IB2I As seen in 
this Figure, this is a system with divided phase space: Sev- 
eral "islands" of closed invariant curves in a stochastic "sea" 
show the existence of an additional conserved quantity other 
than the Jacobian for these orbits. The largest island of 
quasiperiodic orbits in the left half of the surface of section 
corresponds to retrograde orbits, i.e. the stars move around 
the star cluster in the opposite sense to the motion of the star 



